; docformat = 'rst'
;
;+
;
; :Purpose:
;   Aggregate grid cells together into coherant regions based on the method
;   of Rigby et al., (2011) Inversion of long-lived trace gas emissions using 
;   combined Eulerian and Lagrangian chemical transport models, 
;   Atmos. Chem. Phys., 11, 9887-9898, doi:10.5194/acp-11-9887-2011
;
; :Inputs:
;   fp_av: Average footprint
;   q_av: Average emissions
;   lon: Longitudes
;   lat: Latitudes
;   target: Target mole fraction for aggregated footprint*emissions from each region
;
; :Keywords:
;   map: (output, integer) array (dimensions (lon, lat)) of array indices. 0 indicates no emissions.
;
; :Outputs:
;   Hash table of array locations. Table elements are numbered from 0 (no emissions) to N (maximum region index).
;   Each table element is an array of locations.
;
; :Example::
;   IDL> array_locations=mr_aggregate_cells(footprint, emissions, lon, lat, 0.1)
;
;   Output locations of first emissions region:
;   print, array_locations[0]
;     1   2   10   12
;
; :History:
; 	Written by: Matt Rigby and Mark Lunt, University of Bristol, Jul 16, 2013
;
;-
Function mr_aggregate_cells, fp_av, q_av, lon, lat, target, map=map

  compile_opt idl2
  
  nRegions_total=0L

  ;Define some arrays
  ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
  
  LonSize=n_elements(lon)
  LatSize=n_elements(lat)
  
  lon_2d=fltarr(LonSize, LatSize)
  lat_2d=fltarr(LonSize, LatSize)
  
  for LatI=0, LatSize-1 do Lon_2d[*, LatI]=Lon
  for LonI=0, LonSize-1 do Lat_2d[LonI, *]=Lat
  
  area=mr_areagrid(lon, lat)

  ;Find 'center' of footprint
  dummy=max(fp_av, i)
  ai0=array_indices(fp_av, i)
  distarr0=mr_mapdist(lon[ai0[0]], lat[ai0[1]], lon, lat, /nowrap)

  ;Find countries in region
  country=mr_find_country(lon, lat, /ocean)

  ;Calculate average signal that each element contributes to measurement (ppt)
  signal = fp_av*q_av

  ;Array to contain region numbers
  region=long(signal gt 0. and finite(country))
  region[where(region)]=-1L
  
  
  ;Begin creating regions
  ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

  ri=1
  region_total_rec=[0.]
  count=0L
  
  repeat begin
    region_total=0.
    region_total+=max(signal*float(region eq -1), i)
    region[i]=ri
    count++
    ai=array_indices(signal, i)
    
    dist=distarr0[ai[0], ai[1]]
    
    incr=1
    
    while (region_total lt target) do begin

      ar_signal=mr_littleregion(signal, ai[0], ai[1], indices=ai_lr, /nowrap, regionsize=max([9, incr^2]))
      ar_country=mr_littleregion(country, ai[0], ai[1], indices=ai_lr, /nowrap, regionsize=max([9, incr^2]))
      ar_region=mr_littleregion(region, ai[0], ai[1], indices=ai_lr, /nowrap, regionsize=max([9, incr^2]))

      wh=where(ar_region eq -1 and ar_country eq country[ai[0], ai[1]])

      if wh[0] ne -1 then begin
        region_total+=total(ar_signal[wh])
        for n=0L, n_elements(wh)-1 do begin
          region[ai_lr[0, wh[n]], ai_lr[1, wh[n]]]=ri
        endfor
        count++
      endif
      
      if (total(region eq -1) eq 0.) then break
      if incr gt round(dist) then break

      incr+=2

    endwhile

    region_total_rec=[region_total_rec, region_total]
    ri++
  endrep until (total(region eq -1) eq 0.)
  ri_max=ri-1

  ;Check that target is similar order of magnitude
  if abs(alog10(mean(region_total_rec)) - alog10(target)) gt 4 then begin
    print, "Check target order of magnitude: "
    print, " ... target=" + string(target)
    print, " ... mean region size=" + string(mean(region_total_rec))    
    print, " ... returning non-aggregated"
    return, 0
  endif

  ;Tidy up small regions
  ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

  wh=where(region_total_rec lt target/10. and region_total_rec gt 0.)
  while wh[0] ne -1 do begin
    region_total=0.
    wh=wh[sort(region_total_rec[wh])]
    ;Find approx center
    wh_region=where(region eq wh[0])
    ai=array_indices(region, wh_region)
    lonI=round(mean(ai[0, *]))
    LatI=round(mean(ai[1, *]))
    dist=mr_mapdist(LonI, LatI, indgen(LonSize), indgen(latsize), /nowrap)
    ri=1.
    repeat begin
      region_ri=region[where(dist lt ri)]
      ri+=0.5
    endrep until total(region_ri ne wh[0] and region_ri ne 0) gt 0.
    region_ri=region_ri[where(region_ri ne 0 and region_ri ne wh[0])]
    region[wh_region]=region_ri[0]
    region_total_rec[region_ri[0]]+=region_total_rec[wh[0]]
    region_total_rec[wh[0]]=0.
    wh=where(region_total_rec lt target/10. and region_total_rec gt 0.)
  endwhile

  
  ;Re-label regions
  ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
  
  wh=where(region_total_rec gt 0.)
  region_temp=intarr(lonsize, latsize)
  for ri=0, n_elements(wh)-1 do begin
    region_temp[where(region eq wh[ri])]=ri+1
  endfor
  region=region_temp
  ri_max=max(region)

  
  ;Store aggreagted region locations (maybe re-name this to ar_wh?)
  ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
  
  ar_wh=hash()
  for ri=1, ri_max do begin
    wh=where(region eq ri)
    if wh[0] ne -1 then ar_wh[ri-1]=wh
  endfor


  ;Check that emissions still match
  ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
  
  q_agg=0.
  for n=0, n_elements(ar_wh.Keys())-1 do begin
    if total(q_av[ar_wh[n]]) eq 0. then stop
    q_agg+=total(q_av[ar_wh[n]])
  endfor

;  if abs(total(q_av) - q_agg)/q_agg gt 0.01 then begin
;    message, strcompress('Aggregated emissions do not agree agree: ' + $
;      string(total(q_av)) + ' ' + string(q_agg))
;  endif
  
  nRegions_total+=n_elements(ar_wh->keys())

  print, 'CELS_AGGREGATE: Total aggregated regions ', nRegions_total
 
  map=region
 
  return, ar_wh  
  
End
